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1  Report  Summary 

Objective: 

The  objective  is  to  identify  loss  and  degradation  mechanisms  relevant  to  the  frequency  control 
resonator  performance  as  the  resonator  dimensions  reduce  to  the  nanodimensional  spatial  scale. 
The  Department  of  the  Army  interest  in  scaling  such  devices  into  the  nanodimensional  region 
stems  from  the  technical  projection  that  such  resonators  can  operate  at  frequencies  into  the  low 
GHz  spectral  region  thus  providing  a  low  cost,  integratable  NEMS  solid  state  device  option  for 
frequency  control  electronic  applications  relevant  to  jam  resistant  secure  communications.  The 
key  questions  regarding  NEMS  resonator  performance  are  the  fundamental  physical  limitations 
that  arise  in  quality  factor  Q  and  noise  figure  as  the  resonator  is  reduced  to  nanodimensions  in 
order  to  achieve  the  higher  frequencies  of  device  operation.  Such  physical  limitations  may  also 
present  barriers  in  other  NEMS  smart  actuator  and  sensor  electronic  device  applications. 

Approach: 

The  technical  approach  considers  non-equilibrium  heat  generation  and  redistribution  processes 
from  mechanical  strain  during  high  frequency  NEMS  operation  beyond  the  conventional  heat  dif¬ 
fusion  and  local  temperature  approximation.  A  semiclassical  phonon  dynamical  picture  is  intro¬ 
duced  to  go  beyond  the  conventional  models.  Scaling  laws  relevant  to  the  appropriate  phonon 
transport  regimes  and  their  transition  boundaries  are  delineated,  analyzed,  and  compared  with 
results  of  detailed  microscopic  descriptions.  Thus,  the  research  approach  involves  scaling  anal¬ 
ysis,  coupled  with  analytical  and  numerical  modeling  of  phonon  flow  and  heat  redistribution  in 
nanoscale  resonators.  Specifically,  the  following  technical  milestones  were  pursued: 

1 .  Establish  a  picture  of  the  beam  dynamics  in  a  non-equilibrium  thermodynamical  framework. 
This  requires  a  more  refined  look  at  the  approximations  of  the  Euler-Bernoulli  equation 
both  with  and  without  the  Timoshenko  rotational  inertia  and  shear  corrections  to  the  beam 
dynamics,  and  generalizations  of  stress-strain  relations  in  the  presence  of  thermal  gradients. 

2.  For  intrinsic  phonon-mediated  dissipation  mechanisms,  formulate  scaling  laws  based  on  re¬ 
lations  between  the  key  spatial  and  temporal  characteristics  of  phonons,  i.e.,  mean  free  path 
lph  and  phonon  relaxation  time  rph,  and  the  relevant  physical  parameters  of  the  NEMS  struc¬ 
ture,  that  is  beam  thickness  t,  length  L,  and  resonator  operational  frequency  v. 

3.  Scaling  laws  relevant  to  the  appropriate  phonon  transport  regimes  and  their  transition  bound¬ 
aries  will  be  analyzed  and  compared  with  results  of  more  microscopic  analysis  and  theory. 

4.  Analyze  both  surface-assisted  and  anharmonic  coupling  of  flexural  modes  with  thermalized 
phonons,  focusing  especially  on  treatment  of  the  non-equilibrium  phonons  in  the  limit  of  the 
high  frequency  fundamental  mode  of  the  resonator. 
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5.  Formulate  robust  treatment  of  phonons  under  non-equilibrium  thermodynamic  conditions  in 
flexural  nanobeams  by  introducing  Boltzmann  (or  quantum  Wigner-Boltzmann)  framework 
for  description  of  phonon  transport;  conduct  analysis  of  internal  thermal  phonon  flow  inside 
the  beam. 

6.  Consider  how  influences  of  the  naturally  existing  beam  surface  strain  and  surface  roughness 
can  be  integrated  into  the  general  description  of  the  phonon  flow.  Evaluate  surface  effects  and 
their  scaling  properties  in  terms  of  surface-to-volume  importance  as  the  resonator  geometry 
approaches  the  nanolimit. 

7.  Develop  a  reliable  quantification  of  the  Q  factor  and  noise  characteristics  in  presence  of 
non-diffusive  phonon  transport. 

Relevance  to  Army: 

NEMS  have  recently  attracted  great  interest  in  miniaturized  device  design  and  actuation  for  unique 
potential  in  military  and  dual  use  civil  applications.  Their  impact  is  pervasive,  spanning  the  broad 
fields  of  metrology,  engineering,  life  science,  and  medicine.  In  particular,  for  Department  of  Army, 
NEMS  resonators  have  been  the  focus  of  great  attention  for  light-weight,  robust,  on-chip,  high-0 
resonators  and  oscillators  for  frequency  control  electronics  relevant  to  secure  communications. 
Also,  the  possibility  of  obtaining  resonators  with  high  operating  frequencies  into  the  1-10  GHz 
region,  with  small  mass  and  high-0  suggests  NEMS  resonators  and  cantilevers  as  candidates  for 
actuators  and  sensors  in  smart  and  autonomous  systems  with  unprecedented  sensitivities. 

Accomplishments : 

•  Material  parameters  were  successfully  identified  for  the  basic  high  frequency  NEMS  res¬ 
onator  designs;  “operational  frequency  vs.  spatial  dimensions”  maps  have  been  established 
for  the  NEMS  resonator  from  the  Euler-Bernoulli  equation  including  the  Timoshenko  rota¬ 
tional  inertia  and  shear  corrections. 

•  Major  intrinsic  dissipative  mechanisms  related  to  thermoelastic  loss  and  phonon-phonon 
interactions  and  their  scaling  properties  have  been  identified  and  strategically  mapped  to 
provide  insightful  physical  analysis. 

•  From  the  developed  scaling  studies,  it  is  theoretically  noted  that,  in  the  1-10  GHz  operational 
frequency  regime,  the  NEMS  resonator  thermal  dynamics  routinely  goes  beyond  the  limits 
of  the  local  temperature  approximation,  not  at  all  due  to  the  often  considered  time  constraint 
of  the  high  frequency  limit  (which  requires  v  «  r~^,  with  rp/,  «  10  ps  at  room  temperature), 
but,  in  contrast,  due  to  sharp  spatial  inhomogeneity  in  strain  pattern  induced  by  flexure  across 
the  thin  (t  <  lph ,  with  lph  ~  50  nm  at  T  =  300  K  in  Si)  beam  cross  section.  This  spatial 
consideration  leads  to  rapid  ballistic  transfer  of  phonons  across  the  beam  and  suppression 
of  the  dissipation  mechanism  associated  with  the  entropy  production  due  to  inter-branch 
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phonon-phonon  thermal  equilibration.  The  theoretical  analysis  is  formulated  and  conducted 
in  the  Boltzmann  framework  to  capture  the  proper  phonon  dynamics. 

•  Surface  properties  are  heuristically  included  into  the  model  description  opening  avenues  to 
meaningful  analytical  assessment  of  the  surface  effects  on  the  dissipation  in  NEMS  and 
their  effects  on  surface-to-volume  scaling  as  the  resonator  scales  toward  the  ultimate  NEMS 
limit;  as  the  resonator  shrinks,  the  surface  effects  become  influential  in  degrading  the  Q  of 
the  lowest  modes  of  the  resonator. 

Collaborations  and  Technology  Transfer 

•  The  results  of  this  investigation  are  shared  and  discussed  on  the  regular  basis  with  Dr.  Madan 
Dubey  of  ARL  and  Dr.  Arthur  Ballato  of  CECOM. 

•  We  are  in  contact  with  the  group  of  Prof.  M.  Zikry  at  NCSU  for  the  possibility  of  advanced 
multiscale  numerical  modeling  of  these  systems. 

Resulting  Journal  Publications 

A  scientific  publication  describing  heat  dissipation  in  nanoscaled  flexural  resonators  is  anticipated. 
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2  Approach,  rationale,  and  accomplishments 


Below  we  describe  details  of  our  results  and  primary  approach.  The  technical  approach  considered 
non-equilibrium  heat  generation  and  redistribution  processes  from  mechanical  strain  during  high 
frequency  NEMS  operation  beyond  the  conventional  heat  diffusion  and  local  temperature  approx¬ 
imation.  A  semiclassical  phonon  dynamical  picture  was  introduced  to  go  beyond  the  conventional 
models.  Scaling  laws  relevant  to  the  appropriate  phonon  transport  regimes  and  their  transition 
boundaries  were  delineated,  analyzed,  and  compared  with  results  of  detailed  microscopic  descrip¬ 
tions.  Thus,  the  research  approach  involved  scaling  analysis,  coupled  with  analytical  and  numerical 
modeling  of  phonon  flow  and  heat  redistribution  in  nanoscale  resonators. 


2.1  Picture  of  the  beam  dynamics. 

Flexural  modes  available  to  the  double  clamped  resonator  for  large  length  to  thickness  aspect  ratios 
L/t  >  10  can  be  described  analytically  by  the  Euler-Bernoulli  equation  of  solid  mechanics  in  the 
limit  of  no  energy  dissipation,  and  in  the  continuum  approximation  [3,  1].  Applying  this  long- 
beam  approximation  of  the  Euler-Bernoulli  theory,  one  arrives  at  the  equation  for  the  transverse 
displacement  u(z,  t )  for  an  undriven  beam  along  the  direction  of  its  length,  z,  as 


B2u  d 2  ^  B2u  n 

pAW  +  d^EId^~  ' 

Here  p  is  the  material  density,  A  is  the  area  of  the  beam  cross-section,  E  is  the  Young  modulus, 
/  is  the  cross-sectional  area  moment  of  inertia  (bending  moment  of  inertia),  and  El  is  known  as 
the  flexural  rigidity  of  the  beam.  The  boundary  conditions  specify  the  beam  clamping  conditions. 
In  particular,  for  the  rigidly  clamped  boundary  condition,  u  =  0  and  u'  =  du/dz  =  0.  Thus,  for 
the  resonator,  the  eigensolution  for  a  homogeneous  beam  of  length  L  and  rectangular  cross  section 
w  x  t  (i.  e.,  A  =  wt,  I  =  wt3/ 12),  clamped  at  beam  ends  at  z  =  0  and  z  =  L,  can  be  written  in 
phasor  form  as 

un(z,t)  =  Un(z)elUnt , 


where 


Un(z)  =  Cnl[cos(knz)  -  cosh(/cn,z)]  +  Cn2[sm(knz)  -  sinh (knz)]^ 


with  frequency 


I  El. 


kin  =  \  I  —rK  =  \  I  —tk„ 
pA 


(1) 


Here  kn  is  chosen  to  satisfy  the  boundary  conditions,  and  is  a  solution  of  the  transcendental  equa¬ 
tion 

cos (knL)  cosh (knL)  =  1; 
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it  is  found  that  knL  ~  4.730,  7.853,  10.996  for  the  first  three  modes.  Normalization  for  the  mode 
amplitudes  is  conventionally  chosen  as 


L 

dz  U*Um  =  L38nm. 

Similar  equations  can  be  derived  for  the  cantilever  structure,  the  difference  is  in  the  boundary 
condition  at  one  end,  which,  unlike  double  clamped  bridge,  is  left  free  for  the  cantilever. 

While  the  Euler-Bernoulli  theory  works  well  for  the  relatively  long  beams,  the  importance  of 
the  Timoshenko  rotational  inertia  and  shear  corrections  [4]  to  the  beam  dynamics  progressively 
increases  for  smaller  L/t  ratios.  Of  course,  our  consideration  here  is  limited  in  both  approaches  to 
the  small  flexural  vibrations,  before  nonlinearities  start  to  build  up  in  the  resonator  behavior. 

As  a  result  of  our  efforts,  material  parameters  were  successfully  identified  for  the  basic  high  fre¬ 
quency  NEMS  resonator  designs;  “operational  frequency  vs.  spatial  dimensions”  maps  have  been 
established  for  the  NEMS  resonator  from  the  Euler-Bernoulli  equation  including  the  Timoshenko 
corrections. 

In  Fig.  1  we  present  such  a  diagram  showing  frequency  of  the  Si  double  clamped  NEMS  res¬ 
onator  in  dependence  on  the  beam  length  and  thickness.  Two  solid  lines  present  results  of  calcu¬ 
lation  of  the  equi-frequency  fundamental  mode  lines  in  the  framework  of  the  Timoshenko  theory 
allowing  arbitrary  L/t  ratio,  while  dashed  curve  is  a  result  of  the  simple  Euler-Bernoulli  approx¬ 
imation.  Since  the  mode  frequency  in  the  Euler-Bernoulli  framework  is  given  by  Eq.  (1),  the 
dashed  curve  is  simply  a  straight  line  in  the  log-log  plot.  Overall,  for  both  treatments,  the  fun- 


Resonator  length  L  (nm) 


Figure  1:  Graphic  of  NEMS  resonator  operational  frequency  v  =  1  GHz  and  v  =  10  GHz  with  beam 
parameters  length  (L)  and  thickness  (t)  from  the  modified  Euler-Bernoulli  theory. 
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damental  frequency  of  the  oscillator  rises  as  the  resonator  size  and/or  aspect  ratio  (length  to 
thickness)  is  reduced,  as  desired.  In  this  graph  we  have  also  shown  the  phonon  mean  free  path 
(Jph  50  nm  in  Si  at  room  temperature)  on  the  t  axis.  The  following  observations  can  be  made: 

•  Aspect  ratio,  L/t ,  should  be  large  (ideally,  greater  than  10).  Otherwise,  deviations  from 
the  lowest  order  Euler-Bemoulli  theory  increase  rapidly  (e.g.,  via  Timoshenko  corrections); 
this  affects  the  analytics  of  NEMS  resonator:  operational  frequency,  losses,  and  limits  on 
magnitude  of  linear  regime. 

•  10  GHz  resonator  with  a  reasonably  large  aspect  ratio  L/t  (ideally  greater  than  10  or  at 
least  greater  than  5)  will  inevitably  have  characteristically  small  absolute  L  and  t  values. 
In  particular,  its  thickness  can  routinely  be  smaller  than  the  phonon  mean  free  path,  lph, 
establishing  regime  of  ballistic  heat  transfer  in  the  beam  [5].  With  progressive  scaling  of  the 
surface-to-volume  ratio,  surface  effects  will  also  play  a  significant  role. 

•  Local  compression  and  dilation  of  the  beam  material  is  accompanied  by  the  creation  of  non¬ 
equilibrium  phonon  distributions  and  temperature  gradients  and,  consequently,  processes  of 
heat  redistribution  and  entropy  generation.  For  beams  with  substantially  large  aspect  ratios 
(L/t  >  10),  thermal  transport  is  approximately  one-dimensional  (and  takes  place  in  the  plane 
of  the  beam  cross  section);  otherwise,  the  picture  of  heat  redistribution  from  mechanical 
strain  is  complicated  and  multidimensional. 

2.1.1  Structure  of  the  strain  field 

The  structure  of  the  strain  field  corresponding  to  a  particular  flexural  mode  can  be  established  as 
follows.  For  thin  beams,  the  only  nonzero  components  of  the  strain  tensor  are  defined  by  the  mode 
deflection  amplitude  u(z)  as 


d2u 

(2) 

Uzz~  V  dz2' 

^xx  ^ yy  (■ 

(3) 

where  a  is  the  Poisson  ratio.  The  trace  of  the  strain  tensor  is  then 

^Uii  =  ~(1  -  2 <r)yu". 

(4) 

i 


The  strain  pattern  is  visualized  in  Fig.  2  for  the  fundamental  flexural  mode  of  the  double  clamped 
resonator.  An  important  consideration  is  that,  assuming  that  linear  thermal  expansion  coefficient 
of  the  material  is  not  zero,  local  compression  and  dilation  will  be  accompanied  by  the  creation  of 
non-equilibrium  phonon  distributions,  temperature  gradients  and,  consequently  processes  of  heat 
redistribution  and  entropy  generation.  Another  important  observation  is  that  only  for  beams  with 
large  aspect  ratios  (e.  g.,  L/t  greater  than  10),  thermal  transport  is  approximately  one-dimensional. 
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2.1.2  Mode  energy  density  and  total  mechanical  energy 

The  linear  energy  density  along  the  beam  is  given,  in  the  phasor  form,  by 

hw(z)  =  ^  [U(z)f  +  [U”(z)f , 

were  the  first  term  is  the  kinetic  energy  and  the  second  term  is  the  potential  energy  of  flexural 
deformation  of  the  beam.  In  Fig.  3  the  diagram,  the  kinetic  energy  density,  in  phasor  form,  is 
proportional  to  the  square  of  the  fundamental  mode  eigenfunction  Uf  (red  line)  and  the  flexural 
deformation  energy  is  proportional  to  ( U ")2  (green  line).  Both  energy  distributions  are  highly 
uneven  along  the  beam  and  scale  in  accordance  with  L.  However,  as  L  is  reduced,  the  relative 
amounts  of  kinetic  and  flexural  energy  along  the  beam  remain  unchanged  for  the  Euler-Bemoulli 
normalized  eigenmode  U] ;  this  kinetic-flexural  energy  partition  could  vary  with  L  only  by  consid¬ 
ering  higher  order  corrections  to  the  Euler-Bemoulli  equation. 

Finally,  the  total  energy  of  the  mechanical  oscillation  is  given  by 

W  =  J  dzhw(z).  (5) 


2.2  Intrinsic  dissipation  mechanisms  and  their  scaling 

Numerous  mechanisms  were  previously  proposed  as  possible  limitations  to  the  quality  factors  of 
micro-  and  nanoresonators  [6,  7,  8,  9].  Losses  are  of  primary  importance  especially  for  frequency 
control  devices,  as  they  define  noise  characteristics  by  outcome  of  th e  fluctuation-dissipation  the¬ 
orem.  Some  mechanisms  can  be  identified  as  extrinsic  to  the  microresonator  and,  thus,  in  principle 
can  be  eliminated,  while,  intrinsic  mechanisms,  persist  due  to  core  material  properies  of  the  beam. 
For  example,  air  damping  can  be  eliminated  in  the  vacuum-operated  resonators.  Doping  impuri¬ 
ties,  effect  of  the  attachment/support  structure,  and  surface  absorption/desorption  processes  can  be 
accounted  for  as  external. 

For  intrinsic  phonon-mediated  dissipation  mechanisms,  scaling  laws  can  be  formulated  based 


Figure  2:  Time-alterating  compressive  (red)  and  tensile  (blue)  strains  in  the  beam  are  shown  at  deflection 
down  for  the  fundamental  flexural  mode  of  the  double  clamped  resonator.  To  recall,  we  are  interested  in  the 
limit  L/t^$>  1,  L  ss  ultrasubmicron. 
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on  relations  between  the  key  spatial  and  temporal  characteristics  of  phonons,  i.e.,  mean  free  path 
lpi i  and  phonon  relaxation  time  rph,  and  the  relevant  physical  parameters  of  the  NEMS  structure, 
that  is  beam  thickness  t,  length  L,  and  resonator  operational  frequency  v. 

We  start  with  a  short  review  of  the  relevant  phenomena  involved  in  the  heat  transfer  processes 
on  micro-  and  nanoscale. 


2.2.1  Fourier’s  law 

Standard  formulation  of  the  equation  for  the  macroscopic  heat  flow 

q  =  —  kVT 


is  named  Fourier’s  law  and  connects  heat  flow  with  the  temperature  gradient.  The  coefficient  of  the 
proportionality  is  the  thermal  conductivity  k.  The  failure  of  this  formulation  on  the  microscopic 
time  and  space  scales  is  due  to  neglectance  of  the  actual  physical  processes  responsible  for  the  heat 
transfer.  Inclusion  of  these  mechanisms  leads  on  the  macroscopic  level  to  the  appearance  of  the 
microscopic  build-up  time  constant  r, 


q  +  t 


dq 

dt 


-nWT 


the  formulation  typically  attributed  the  Cattanco  and  Vernotte  who  independently  proposed  this 
heuristic  generalization.  This  leads  to  the  so-called  hyperbolic  heat  equation 


”  dt  v  dt 2 


=  —  p  —  T 


dp 
dt ' 


Figure  3 :  Linear  energy  density  in  phasor  form  for  the  fundamental  flexural  mode  of  the  double  clamped 
resonator. 
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Here  Cv  is  the  specific  heat,  p  is  the  power  density.  In  absence  of  r  we  obtain  the  classical  heat 
diffusion  law.  Explicitly, 

3  K 

T  =  c^2 

is  the  microscopic  thermal  time.  On  the  macroscopic  level,  it  brings  up  the  phenomenon  of  the 
non-diffusive  heat  transfer  by  means  of  the  quickly  decaying  heat  waves. 

Taking  for  silicon  at  room  temperature  the  standard  set  of  parameters  c  =  5880  m/s=  5.88  nm/ps, 
k  =  160  W/mK,  and  Cv  =  1.78  x  105  J/m3K,  we  obtain  for  the  characteristic  phonon  microscopic 
time  tpi i  ~  3x3.5  ps~  10  ps,  which  corresponds  to  the  phonon  mean  free  path  lph  =  cTph  ~  50  nm. 

A  different,  macroscopic  time  can  be  associated  in  NEMS  resonator  with  the  heat  diffusion 
across  the  beam  width  as  a  result  of  the  flexural  vibrations.  A  simple  treatment  due  to  Zener  [8] 
reveals 

_  t2Cv  _  3  ft\2  1 
Tdlff  n2n  n2\c)  rph 

Introducing  Tbaiustic  =  \/3t/c  —  a  time  neccesary  for  a  typical  phonon  to  ballistically  traverse  the 
beam  width,  we  conclude  that 


Tdiff 


r2 

ballistic 


7T2  Tph 


3  t 2 


7 T2  l 


2  TPh ' 
ph 


Thus,  quite  obviously,  for  t  A>  lph  slow  diffusive  heat  relaxation  dominates  and  the  heat  transfer  is 
essentially  Fourier’s  transfer. 

The  standard  theory  of  thermoelastisity  relies  essentially  on  the  validity  of  the  Fourier  heat 
law  [8,  16],  thus,  either  this  theory  should  be  upgraded  to  incorporate  non-Fourier’s  effects  in  the 
calculational  framework  (the  path  chosen,  in  particular  by  Sun  et  al.[17];  it  is  worth  to  mention 
that  estimates  suggest  that  in  the  interval  of  parameter  values  considered  in  this  paper  an  effect  of 
non-Fourrier  heat  transfer  is  minimal),  or  microscopic  relaxation  processes  should  be  considered 
and  accoounted  for  separately  in  the  analysis  of  intrinsic  losses.  In  case  of  insulators  at  elevated 
temperatures,  the  specific  type  of  phonon-phonon  interactions  represent  themselves  as  a  so-called 
Akhiezer  effect. 


2.2.2  Akhiezer  effect 

Akhiezer  effect  is  most  commonly  considered  in  relation  to  the  losses  of  ultrasonic  waves  in  di¬ 
electric  crystals  at  elevated  temperatures  and  is  a  result  of  the  phonon-phonon  interactions.  Due  to 
anharmonicity,  strain  modulates  phonon  frequencies  (and,  strictly  speaking,  in  a  different  fashion 
for  different  phonon  modes).  As  a  result,  original  local  equilibrium  phonon  distribution  becomes 
distorted  and  requires  a  microscopic  time  r  to  reestablish  phonon  equilibrium  locally.  This  irre¬ 
versible  process  requires  absorption  of  elastic  energy  to  generate  entropy  in  the  phonon  sybsys- 
tem.  A  simple  qualitative  treatment  was  published  by  Bommel  and  Dransfeld  [13].  Essentially, 
they  have  split  phonon  modes  in  two  groups  -  those,  experiencing  positive  temperature  changes 
and  those,  experiencing  negative  or  no  changes.  At  the  same  time,  no  intragroup  deviation  from 


10 


the  equilibrium  distribution  was  assumed.  Thus,  intergroup  temperature  relaxation  was  the  main 
process,  giving 


Q'1 


CvT'f2 

E 


fir 


at  sufficiently  low  frequencies  of  the  acoustic  vibrations  fir  <C  1.  Here  7  =<  7^  >  is  the  average 
Griineisen  constant 

(A  c/c)jk  =  7 jk{&p/p), 


defined  by  the  material  anharmonicity  (as  given  by  third-order  elestic  moduli,  or,  e.  g.,  related  to 
the  temperature  expansion  coefficient  aT  as  7  =  3aTK/Cv). 

It  is  interesting  to  see  that  expressing  7  via  a-p  and  K  via  E  ~  3 K  in  the  prefactor  for  micro¬ 
scopic  losses  and  making  use  of  Cv  m  Cp,  we  readily  recover  the  conventional  thermoelastic 
prefactor  ajTE/Cp.  Of  course,  a  number  of  approximations  led  to  this  coinsidence,  in  particular, 
assumption  of  the  high  temperature  T  >  0£>,  which  is  obviously  questionable  at  room  tempera¬ 
ture  for  the  silicon  with  (-)  D  =  665  K,  assumption  of  the  simple  Debye  phonon  spectrum,  and  the 
validity  of  the  application  of  the  Griineisen  constant  average. 

More  advanced  microscopic  treatment  was  proposed  by  Woodruff  and  Ehrenreich  [12],  we 
return  to  this  publication  later. 


2.2.3  Scaling  of  losses  with  resonator  dimensions 

Applying  results  of  these  early  developments,  we  have  analyzed  major  intrinsic  dissipative  mech¬ 
anisms  related  to  thermoelastic  loss  [10]  and  phonon-phonon  interactions  [11].  Their  scaling 
properties  have  been  identified  and  strategically  mapped  to  provide  insightful  physical  analysis. 
Specifically,  from  the  developed  scaling  studies,  it  was  theoretically  noted  that,  in  the  1-10  GHz 
operational  frequency  regime,  the  NEMS  resonator  thermal  dynamics  routinely  goes  beyond  the 
limits  of  the  local  temperature  approximation,  not  at  all  due  to  the  often  considered  time  constraint 
of  the  high  frequency  limit  (which  requires  v  «  t/jP  ,  with  rph  ~  10  ps  at  room  temperature),  but, 
in  contrast,  due  to  sharp  spatial  inhomogeneity  in  strain  pattern  induced  by  flexure  across  the  thin 
(t,  <  lph,  with  lph  ~  50  nm  at  T  =  300  K  in  Si,  see  Fig.  4)  beam  cross  section.  This  spatial 
consideration  leads,  qualitatively,  to  rapid  ballistic  transfer  of  phonons  across  the  beam  and  sup¬ 
pression  of  the  dissipation  mechanism  associated  with  the  entropy  production  due  to  inter-branch 
phonon-phonon  thermal  equilibration.  It  is  further  noted  that  the  more  consistent  theoretical  analy¬ 
sis  should  be  formulated  and  conducted  in  the  Boltzmann  framework  to  capture  the  proper  phonon 
dynamics. 


2.3  Phonons  under  non-equilibrium  thermodynamic  conditions  in  flexural 
nanobeams 

A  quantitative  description  of  the  processes  of  heat  generation  and  redistribution  in  thin  flexural 
beams  has  been  achieved  by  applying  the  Boltzmann  formalism  to  capture  dynamics  of  the  non¬ 
equilibrium  phonon  distribution  resulting  from  the  spatially  inhomogeneous  mechanical  compres- 
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sion  and  dilation.  The  multidimensional  integro-differential  equation  has  been  derived  for  the 
phonon  distribution  function,  encompassing  effects  of  the  modulation  of  the  phonon  frequencies 
with  applied  spacially  inhomogeneous  strain,  phonon  ballistic  transfer,  and  scattering  processes, 
leading  to  thermalization  in  the  phonon  ensemble.  This  allowed  a  reliable  quantification  of  the 
resonator  quality  factor  Q  in  presence  of  non-diffusive  phonon  transport. 


2.3.1  Phonon  Hamiltonian  in  presence  of  strain 

The  phonon  Hamiltonian  for  a  particle  with  wave  vector  q  (for  simplicity  of  notation,  q  includes 
both  the  wave  vector  and  the  index  of  the  phonon  branch)  is 

H  =  hw(q,  r,  t)  =  [1  +  a(q,  r,  t)]  . 

In  the  presence  of  strain,  in  the  lowest  order  [11] 

«(q,r  ,t)  =  -^yik(q)uik(r,t), 

ik 

where  u ik  is  the  strain  tensor  and  (q)  is  the  generalized  Griieneisen  tensor  [12].  The  only  non¬ 
zero  components  of  the  strain  tensor  in  the  beam  were  given  in  Eqs.  (2),  (3)  above.  With  that,  we 
can  write  explicitly 


a(q,r,t) 


(1  +  <7)7*a(q)  -  ^ 

i 


yu"{z). 


Now  identifying  all  diagonal  components  7 a  of  the  tensor  with  the  single  GrUneisen  constant  7  in 
the  theory  of  the  thermal  expansion ,  we  can  simplify  further,  so  that 

a(q, r,  t)  =  -tXI  Uii  =  t1  _  2crhyu"(z)-  (6) 


This  simplification  that  7  is  scalar  and  ^-independent  was  historically  made  in  the  early  devel¬ 
opment  of  bulk  thermal  expansion  theory  and  has  served  fairly  well  in  simple  situations  requiring 
exclusively  phonon  bulk  averages;  in  the  limiting  case  where  the  nanoregime  is  in  effect,  averaging 
with  respect  to  JA  and  yzz  is  not  so  easily  justifiable,  especially  when  one  needs  to  accomo¬ 
date  the  physically  dictated  inherent  anisotropy  into  the  analysis  of  the  degradation  processes.  For 
example,  reducing  to  a  single  constant,  7,  artificially  eradicates  the  Akhiezer  effect.  Therefore,  to 
preserve  the  physical  consequences  of  the  Akhiezer  phenomenology,  especially  at  the  NEMS  spa¬ 
tial  and  frequency  scale,  where  the  effect  is  particularly  important,  one  should  allow  for  at  least  two 
different  constants  71  and  72  to  mimic  the  variation  in  tensor  7^(q).  A  simple  qualitative  treatment 
of  the  Akhiezer  effect  is  then  possible,  as  originally  provided  by  Bommel  and  Dransfeld  [13].  As 
their  starting  point,  they  have  split  all  phonon  modes  in  two  distinguished  groups  —  those,  experi¬ 
encing  positive  temperature  changes  and  those,  experiencing  negative  or  no  temperature  changes. 
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2.3.2  Boltzmann  equation 


Our  quantitative  description  of  the  processes  of  heat  generation  and  redistribution  in  thin  flexu¬ 
ral  beams  has  been  achieved  by  applying  the  Boltzmann  formalism  to  capture  dynamics  of  the 
non-equilibrium  phonon  distribution  resulting  from  the  spatially  inhomogeneous  mechanical  com¬ 
pression  and  dilation.  Introducing  the  instant  distribution  function  N (q,  r,  t)  for  the  population  of 
the  phonons  in  the  beam,  a  thermal  equlibrium  distribution  function 

N0(uj,T)  =  N0(hiu/kT )  =  [exp (hu/kT)  -  l]-1 , 


its  derivative  N{ q  in  respect  to  x  =  hu/kT,  further  splitting  N  into  the  “equilibrium”  part  N0(hu/kT0), 
corresponding  to  the  local  phonon  frequencies  u  (defined  by  the  spatially  inhomogeneous  strain) 
and  the  equlibrium  temperature  T0,  and  the  non-equilibrium  part  Ar,  (q.  r,  t)  —  N (q,  r,  t)—N0(huj/kT0), 
we  can  linearize  the  original  Boltzmann  equation  for  the  phonon  distribution  function 

/  dN\  d  N  dN  di o  dN  du 

V  )  coil  (Jr  dq  dq  dr 


to  arrive  at 


dN  \  dN\  ,  huj0  da  dN\ 

ar- 


(7) 


Here  the  mode  velocity  v0  =  dco0/dq.  We  would  like  to  note  that  this  particular  linearization  is 
different  from  another  commonly  performed  splitting  N  =  N0(u0,  T0)  +  Nu  where  frequency  u0 
is  used  in  the  true  equilibrium  part. 

Treating  the  collision  term  in  the  effective  relaxation  time  approximation,  we  let 


(dN 

\dt 


coll 


N-N0{u),T) 

'Tph 


(8) 


where  rph( q)  is  the  phonon  relaxation  time.  There  are  a  number  of  subtle  issues  associated  with 
this  equation.  First,  it  is  a  heuristic  approximation  to  the  realistic  microscopic  description  for  the 
phonon-phonon  interactions  [14].  Even  then,  actually  two  classes  of  phonon-phonon  scattering 
processes  take  place,  normal  processes  that  conserve  total  momentum  of  the  colliding  pair  and 
umklapp  process,  that  do  not  conserve  total  momentum.  As  written  in  Eq.  (8),  it  is  implicitly 
assumed  that  the  umklapp  processes  dominate,  so  as  the  asymptotic  phonon  distribution  can  be 
characterized  by  a  single  parameter  —  local  equlibrium  temperature  T.  Further,  this  local  equlib¬ 
rium  temperature  T  should  not  be  defined  by  averaging  phonon  distribution  only  at  a  single  spatial 
point;  instead,  it  should  be  obtained  by  averaging  over  the  area  of  the  size  of  about  lph.  Thus,  it  is 
clear,  that  T  should  approach  T0  as  the  beam  thickness  becomes  smaller  than  lph.  Nevertheless,  it 
can  be  shown  that  treating  T  as  a  single-point  quantitity  produces  virtually  identical  results  for  the 
losses  even  in  the  limit  of  t  <  lph. 

Now  we  consider  a  case  of  substantially  different  spatial  scales  L  3>  t  for  the  Euler-Bernoulli- 
Boltzmann  resonator,  allowing  to  assume  that  all  thermal  transfer  is  essentially  one-dimensional 
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and  takes  place  along  axis  y  in  the  plane  of  the  cross  section.  Seeking  oscillatory  solutions  of  TV) 
with  flexural  fundamental  frequency  fl,  we  obtain 


(1  -  iVtTphjNx  -  v0Tph  cos  9N[ 


(ittTpha(y) 


(9) 


with  the  temperature  difference  T) 
processes,  i.  e., 


=  T  —  T0  defined  by  the  total  energy  balance  in  scattering 


Further  progress  can  be  made  (i)  by  assuming  the  Debye  model  (u  =  c0q )  with  the  ^-independent 
velocity  v0  =  c0  and  (ii)  by  assuming  that  the  relaxation  time  rph  is  also  phonon  mode-independent. 
Expanding  then  Eq.  (10)  to  the  first  order  in  the  relaxation  approximation,  we  obtain 


[  d9  sin  9  [  dq  q2ftw(q)Ni  =  —  [  dqq2 

4:71  /  27T  / 


fa^o(g) 

kTo 


kN'. 


(11) 


The  equation  for  r.h.s.  is  easily  recognizable:  recalling  that  the  specific  heat  of  mode  q  is 


S(q)  = 


kujp (q) 
kTn 


kN', 


and  the  heat  capacity  of  the  phonon  subsystem  (per  unit  volume) 


C  = 


1 

2tt 


I 


dq  q2S(q), 


we  simply  get  for  the  r.h.s.  of  Eq.  (11),  CT^. 

Now  assuming  momentarily  that  the  Griineisen  tensor  7  is  also  a  ("/-independent  scalar  constant 
as  given  by  Eq.  (6),  we  can  multiply  Eq.  (9)  by  phonon  energy  and,  integrating  it  over  modulus  q, 
get  an  equation  for 

f(y,cosd)  =  j- 

(1  -  iflTph)f  -  lph  cos  Of  =<N[T1-  iriTphT0a(y)} .  (12) 

Using  Eq.  (6)  expressing  a(y)  via  7  and  the  flexural  mode  deflection  amlitude  u(z),  it  is  also 
convevient  to  introduce  a  function  /7: 


j  dq  q2huoNi 


f  =  (1  -  2(t)7«"(-)/7- 


Here  the  advantage  lies  in  the  explicit  factorization  of  the  z  dependence,  i.  e.  /7  is  functionally  the 
same  along  the  beam  length.  The  Eq.  (12)  should  be  accompanied  by  the  physically  meaningful 
boundary  conditions  and  then  solved  numerically  either  directly  or  as  a  series  expansion.  To  get  a 
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flavor  of  the  involved  mathematics,  we  mention  a  particular  possible  expantion  of  / 

/  =  X)  Yn{y)Pn{ COS0)  =  r0  +  Vi  COS0  +  .  .  . 

n 


where  Pn(x)  are  Legendre  polynomials.  Making  use  of  their  properties,  one  can  multiply  Eq.  (12) 
for  /  by  various  Pm  and  integrate  over  angle  9  to  obtain 

(7i  77/  1  \ 

+  2^3yn+i  )  =  (i^rphCT0a(y)  +  Y0)  Sn0, 

which  should  be  furhter  solved  accompanied  by  appropriate  boundary  conditions. 


2.3.3  Phonon  related  losses  in  flexural  beams 


As  it  was  mentioned  above,  a  conseptually  similar  problem  of  sound  attenuation  in  dielectric  crys¬ 
tals  was  considered  in  [12,  13].  Unlike  these  publications  where  an  infinite  domain  with  solutions 
in  the  form  of  plane  waves  was  analyzed,  here  we  are  forced  to  operate  with  solutions  on  the  finite 
domain  in  presence  of  spacial  inhomogeneities.  The  average  rate  at  which  the  energy  is  removed 
from  the  flexural  vibration  is  equal  in  the  steady  state  to  the  rate  Z  of  transfering  energy  from  the 
phonon  system  to  the  heat  bath  [15] 


Z 


Using  Eq.  (7)  to  express  the  collisional  term  via  time  derivative  of  N,  drift  and  diffusion  terms, 
and  further  removing  terms  that  time  average  to  zero,  we  find 


Z  = 


167 r3 


>v 


dct 

dr  I  dq  ftw0Re  N*  —  . 


(13) 


With  Q~l  defined  as  the  fraction  of  the  oscillation  energy  lost  per  radian  of  flexural  vibration,  total 
flexural  energy  W  given  by  Eq.  (5),  and  average  rate  of  energy  transfer  Z  defined  by  Eq.  (13),  we 
arrive  at 


Q~l  = 


QW 


for  the  inverse  of  the  resonator  quality  factor  Q. 

Invoking  previously  made  approximations  of  (i)  a  large  beam  aspect  ratio,  (ii)  the  Debye  model, 
and  (iii)  constant  scalar  Griineisen  coefficient,  we  can  perform  a  number  of  partial  integrations  to 
get,  in  terms  of  function  /7: 


z  =  n 


(1  -  2a)2ry2w 
47 r2 


dz  [u'\z)f  dyy  /  dO  sin  9  Irn  {/7} 


for  the  particular  beam  mode  u(z).  With  (or  /)  found,  phonon-related  losses  in  the  beam  can 
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be  numerically  evaluated. 

Again,  we  stress  here  that  when  the  assumption  of  the  one  dimensional  thermal  transport  (case 
of  a  relatively  long  beam,  that  is  L  t)  is  not  valid,  A',  cannot  be  factorized  to  effectively 
separate  the  coordinate  z  along  the  beam,  then  leaving  us  with  a  complex  system  of  PDEs  to  be 
solved  numerically. 

2.3.4  Analysis  of  internal  thermal  phonon  flow  inside  the  beam  in  the  (comparatively  sim¬ 
ple)  limit  of  large  aspect  ratio  L/t  >  1  (L  «  ultrasubmicron) 

Major  outcomes  of  the  phenomelogical  descriptions  can  be  reproduced,  and  corresponding  char¬ 
acteristic  parameters  can  be  established  in  the  framework  of  the  unified  formalism: 

•  Fourier’s  heat  transfer  at  long  times  (t  >  rph)  and  large  distances  (Z  >  lth)  [8,  16]; 

•  non-Fourier  heat  transfer  at  short  times  and  distances  [17];  specifically,  since  NEMS  with 
GHz  frequencies  are  very  small  (t  <  lph),  thermoelastic  mechanism  is  dominated  by  the 
non-Fourier  heat  transfer,  but,  in  any  case,  is  of  relatively  minor  importance; 


Resonator  frequency  v  (GHz) 

Figure  4:  Magnitude  of  losses.  Color  code:  blue — high  Q,  low  losses,  red — low  Q ,  high  losses.  Here 
we  assume  equality  of  the  thermoelastic  and  microscopic  phonon-phonon  prefactors,  and  plot  the  quality 
factor  for  the  first  flexural  mode  of  the  NEMS  resonator  as  a  function  of  the  beam  length  L  and  thickness  t 
due  to  both  mechanisms  together.  To  stay  mostly  in  the  limits  of  validity  of  the  Euler-Bernoulli  theory,  one 
should  consider  only  beams  with  L  >  5 1.  While  prefactor  and  rp/(  depend  only  on  the  material  properties, 
Hi  =  '2itv\  and  rfujj  are  obviously  functions  of  the  beam  geometry.  As  one  reduces  beam  dimensions, 
he  can  clearly  see  a  crossover  in  the  scaling  law  around  t  =a  few  lph,  when,  with  increasing  resonator 
frequency  the  thermoelastic  mechanism  gives  in  to  the  phonon-phonon  microscopic  losses.  Quantitatively, 
for  a  specific  beam  aspect  ratio  L/t  =  10,  a  10  GHz  NEMS  has  substantially  higher  losses  than  a  1  GHz 
resonator. 
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•  distortion  by  strain  of  the  inter-branch  local  thermal  equilibrium  and  inter-branch  phonon- 
phonon  thermalization — the  Akhiezer  dissipation  mechanism  associated  with  the  entropy 
production  due  to  inter-branch  phonon-phonon  thermal  equilibration  [11,  13]  actually  dom¬ 
inates  in  small  structures. 

•  In  quantitative  terms,  a  10  GHz  silicone  NEMS  resonator  with  L  =  10/  has  substantially 
higher  losses  than  a  1  GHz  resonator. 


2.4  Surface  effects 

Rather  common  outcome  of  numerical  simulatons  and  experimental  observations  in  MEMS  is  that 
the  reduction  of  the  resonator  size  leads  to  a  substantial  reduction  in  the  quality  factor  Q  [18,  19]. 
A  feasible  explanation  for  this  trend  is  that  microresonator  mechanical  losses  are  determined  to 
a  great  degree  by  micro  structure  surfaces  and  surface-assisted  processes.  With  down-scaling  of 
the  resonator,  a  larger  fraction  of  the  beam  atoms  are  found  at  the  surface.  Indeed,  a^S/V  grows 
linearly  with  reduction  of  the  beam.  A  word  of  caution  should  be  given,  although:  as  it  was  shown 
above,  progressive  miniaturazation  directly  leads  to  the  increase  in  the  frequency  of  the  fundamen¬ 
tal  mode  and,  consequently,  growth  of  the  losses  due  to  the  Akhiezer  effect.  Coincidentally,  the 
scaling  law  for  this  mechanism  also  follows  the  S/V  dependence. 

Experimental  work  by  Mihailovich  and  MacDonald  [18]  suggest  a  presence  of  losses  due  to 
irregularities  and  damage  of  the  crystal  structure  in  the  vicinity  of  the  reactive-ion  etched  surfaces. 
Obviously,  their  relative  influence  will  extend  with  reduction  of  the  resonator.  Damage-healing 
treatment  was  suggested  to  overcome  this  negative  effect,  considering  that  the  NEMS  structure 
can  withstand  them.  It  was  shown  that,  for  example,  thermal  oxidation  can  reduce  losses  due  to 
etched  surfaces  by  a  factor  of  two. 

Surface  tension  and  the  possible  reconstruction  on  the  surface  lead  to  the  size-dependence  of 
the  effective  elastic  modulii  and  material  density,  the  phenomenon  observed  experimentally  [20], 
as  well  as  in  the  numerical  simulation  [19].  Authors  of  [19]  even  suggest  that  they  were  also 
able  to  identify  effects  of  the  reconstruction  of  the  material  at  corners  of  the  beam.  These  effects 
of  surface  on  the  effective  material  constants  directly  contribute  to  the  further  deviation  from  the 
established  scaling  laws  for  the  resonator  frequency,  but  only  indirectly  and  negligibly  to  the  loss 
factors.  A  different  mechanism  —  notable  enhancement  of  the  third-  and  fourth-order  anharmonic 
phonon-phonon  interaction  in  presence  of  the  static  strain  in  the  surface  layer  could  leads  to  the 
direct  dependence  of  the  r^am  on  the  surface-to- volume  ratio.  For  example,  one  can  qualitatively 
write 

1  1  aoS  1 

-  - 1 - 

rr beam  bulk  T/  surf 

Tph  Tph  17  Tph 

where,  presumbly,  Tp}'^  <C  Tp%lk. 
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2.4.1  Boundary  conditions  at  beam  interfaces  for  phonon  distribution  function 

We  have  developed  a  heuristic  approach  to  incorporate  both  surface-assisted  and  anharmonic  cou¬ 
pling  of  flexural  modes  with  thermalized  phonons,  focusing  especially  on  treatment  of  the  non¬ 
equilibrium  phonons  in  the  limit  of  the  high  frequency  fundamental  mode  of  the  resonator.  As 
a  result,  influences  of  the  naturally  existing  beam  surface  strain  and  surface  roughness  could  be 
integrated  into  the  general  description  of  the  phonon  flow  through  a  set  of  heuristic  parameters  de¬ 
scribing  scattering  and  thermalization  due  to  interfaces.  For  example,  one  can  heuristically  force 
the  phonon  distribution  N  to  satisfy  the  physically  meaningful  boundary  condition  for  the  partial 
specular  reflection,  i.  e., 

N\n_e  =  pN |0<7r/2  +  (1  -  p)Nq(u>,  Tb),  (14) 

or,  linearized  in  terms  of  the  N\ 

j-  rjn 

mu-g  =  pn^/2  -  (i  -  P)K—^,  (i5) 

Here  the  parameter  0  <  p  <  1  quantifies  the  part  of  reflected  phonons:  p  is  unity  for  perfect 
specular  reflection  and  zero  for  no  specular  reflection  contribution;  we  further  assume  that  the 
complementary  part,  (1  —  p),  is  instantly  thermalized  due  to  interaction  with  the  interface — this 
assumption  is  actually  slightly  different  than  the  usual  assumption  of  isotropic  diffusive  scattering 
[21].  Further,  Tb  =  T0  +  Tbl  is  the  local  boundary  temperature.  This  approach  allows  for  incor¬ 
poration  of  surface  effects  and  estimation  of  their  scaling  properties  in  terms  of  surface-to-volume 
importance  as  the  resonator  geometry  approaches  the  nanolimit. 

It  is  noted  that  the  influence  of  “total”  thermalization  (1  —  p)  at  interfaces  can  be  qualitatively 
interpreted  in  the  following  way  —  assuming  dominating  ballistic  transfer  of  phonons  across  the 
width  of  the  beam,  then  it  follows  as  a  consequence  that 

1  1  1  -p 

-  ~ - 1 - 

rj-beam  __ bulk  ’ 

•  ph  •  ph,  '  ballistic 

which  means  that  the  position  of  the  phonon  loss  peak  in  the  lower  part  of  Fig.  4  shifts  away 
from  its  original  vertical  location  at  approximately  ^lr^lk  =  1  as  a  function  of  the  specularity 
coefficient  p.  Therefore,  one  can  expect  an  influence  of  p  on  the  intensity  of  degradation  processes 
and,  therefore,  Q  in  the  high  frequency  limit. 

From  this  vignette,  one  can  see  an  excellent  opportunity  to  explore  further  the  overall  impor¬ 
tance  of  surface  connected  behavior  on  scaling  properties  of  NEMS  structures  in  a  user  friendly 
and  quantitative  manner. 

2.5  Final  observations 

In  summary,  we  have  considered  non-equilibrium  heat  generation  and  redistribution  processes 
from  mechanical  strain  during  high  frequency  NEMS  operation  beyond  the  conventional  heat  dif- 
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fusion  and  local  temperature  approximation.  A  semiclassical  phonon  dynamical  picture  was  in¬ 
troduced  to  go  beyond  the  conventional  models.  Scaling  laws  relevant  to  the  appropriate  phonon 
transport  regimes  and  their  transition  boundaries  were  delineated  and  analyzed;  the  advanced  the¬ 
oretical  analysis  was  formulated  and  conducted  in  the  Boltzmann  framework  to  capture  the  proper 
phonon  dynamics  inside  flexural  beams.  Further,  surface  properties  were  heuristically  included 
into  the  model  description  opening  avenues  to  meaningful  analytical  assessment  of  the  surface 
effects  on  the  dissipation  in  NEMS. 

As  a  result  of  this  study,  the  intrinsic  mechanisms  and  related  analytic  principles  defining 
NEMS  degradation  with  scaling  have  been  developed,  and  have  been  expressed  in  a  user  friendly 
form  for  application.  It  is  timely,  therefore,  to  consider  a  follow  on  study  to  pursue  the  analytic  and 
numerical  application  of  the  developed  principles  so  that  realistic  values  of  Q  can  be  calculated  for 
actual  resonator  materials  and  geometries. 
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